Experimentally Constrained Molecular Relaxation: The Case of Glassy GeSe2 
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An ideal atomistic model of a disordered material should contradict no experiments, and should 
also be consistent with accurate force fields (either ab initio or empirical). We make significant 
progress toward jointly satisfying both of these criteria using a hybrid reverse Monte Carlo approach 
in conjunction with approximate first principles molecular dynamics. We illustrate the method by 
studying the complex binary glassy material g-GeSe2. By constraining the model to agree with 
partial structure factors and ab initio simulation, we obtain a 647-atom model in close agreement 
with experiment, including the first sharp diffraction peak in the static structure factor. We compute 
the electronic state densities and compare to photoelectron spectroscopies. The approach is general 
and flexible. 

PACS numbers: 61.43.Fs, 71.23.Cq, 71. 15. Mb, 71. 23. An 
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INTRODUCTION 

The modeling of complex materials based upon molec- 
ular dynamics simulation has been one of the remark- 
able advances in theoretical condensed matter physics. 
Whether the potentials chosen are empirical or ab initio, 
remarkable insights have accrued for diverse problems in 
materials physics and beyond. There is, however, an un- 
satisfying point to the logic of MD simulation: it does not 
make use of all the information available about a material 
under study - notably the information implied by experi- 
ments. Simulations often cannot achieve agreement with 
experiment because of short simulation times, small sys- 
tem sizes or inaccuracies in the interactions. Successful 
prediction of new properties is more likely for models in 
agreement with existing data. Imposition of experimen- 
tal information may be important in phase-separated or 
other complex materials for which obtaining a suitable 
starting structure may be difficult, and for which short 
MD time scales preclude the emergence of such structures 
in the model. 

A different approach to model construction imple- 
mented by McGreevyd S and colleagues is the 
so-called "Reverse Monte Carlo" (RMC) method. Here, 
one explicitly sets out to make an atomistic model which 
agrees with experiments. RMC has been widely used to 
model a variety of complex disordered materials. This 
is accomplished by making Monte Carlo moves which 
drive a structural model toward exact agreement with 
one or more experiments. In practice, RMC is the ideal 
method to explore the range of configurations which are 
consistent with experiment (s). Without adequate limi- 
tation to a "physical" subspace of configuration space, 
it is unlikely to produce a satisfactory model. That is, 
only a subset of RMC models [which match the experi- 
ment(s)] is physically realistic (consistent with accurate 
interatomic interactions). The imposition of topologi- 
cal/chemical bonding constraints in RMC can ameliorate 
this problem, but not remove it entirely (sl. The math- 



ematical structure of constrained RMC is a constrained 
optimization "traveling salesman" problem. In our pre- 
vious implementation of constrained RMC we formed a 
positive definite (quadratic) cost or "penalty" function ^, 
which was then minimized (ideally, but not practically) 
to zero for a structural model which exactly satisfies all 
constraints imposed: 

K Mj L 
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where rf^ is related to the uncertainty associated with 
the experimental data points, K is the number of exper- 
imental data sets employed, Mk is the number of data 
for the K^^ set and L is the number of additional (non- 
experimental) constraints included. The quantity Q is 
the appropriate generalized variable associated with ex- 
perimental data i^(Q) and P; > is the penalty function 
associated with each additional constraint and A/ > 0. 
Such "additional" constraints can be of many different 
forms (for example, one may impose chemical or topo- 
logical ordering, or phase separated units within a con- 
tinuous random network) The coordinates of atoms are 
changed according to Monte Carlo moves, which is akin 
to a simulated annealing minimization of our cost func- 
tion ^. The method is easy to implement, though care 
must be taken to include the minimum number of inde- 
pendent constraints possible to reduce the likelihood of 
getting "stuck" in spurious minima. We have shown that 
inclusion of suitable constraints leads to models of a-Si 
much improved compared to RMC models using oiily the 
structure factor (first term of Eq. ^) as constraint pj ■ 

As the creation of models of complex materials is a dif- 
ficult task, it is of obvious advantage to incorporate all 
possible information in fabricating the model. We assert 
that an ideal model of a complex material should (1) be 
a minimum (metastable or global) of a suitable energy 
functional faithfully reproducing the structural energet- 
ics, (2) should contradict no experiments. When stated 
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in these terms, our criterion seems quite obvious, yet cur- 
rent simulation schemes do not simultaneously accommo- 
date both criteria, but focus only on one or the other. 

In this paper, we merge the power of ab initio molec- 
ular simulation with the a priori information of experi- 
ments to create models consistent with experiments and 
the chemistry implied by accurate interatomic interac- 
tions. To obtain joint agreement, we unite MD with 
the Reverse Monte Carlo (RMC) method. We name 
the scheme "Experimentally Constrained Molecular Re- 
laxation" (ECMR). One can understand our scheme as 
a way to "tune" a structural model using MD within 
the space of experimentally realistic models as defined by 
RMC. We choose a troublesome and complex material 
with a long experimental and modeling history: g-GeSe2. 

From an algorithmic perspective, our scheme has im- 
portant advantages. For example, to model a glass like 
GeSe2 or Si02 using first principles methods, the method 
of choice is to form an equilibrated liquid, use some dis- 
sipative dynamics to simulate an (unphysically) rapid 
quench of the liquid into an arrested phase and finally 
to relax this to a local energy minimum, presumably 
at astronomical Active temperature (high potential en- 
ergy) . Usually some repeated "annealing" cycles are also 
used. If one is interested in a glassy phase all the work 
of forming and equilibrating the liquid is redundant, and 
it is a pious hope that the arrested liquid will resemble 
a real glassy phase. Evidently the likelihood for suc- 
cess is strongly affected by topological and chemical sim- 
ilarity of the melt to the physical amorphous phase. If 
complex ordering "self-organization" or phase separation 
occurs in the physical amorphous phase, the short simu- 
lations of conventional ab initio schemes will surely miss 
these important structural features. In this vein, we have 
used ECMR to construct models of a-Si with intermedi- 
ate range order on a nanometer length scale by in- 
clusion of Fluctuation Electron Microscopy 0- We note 
that successful techniques do exist to tackle the time- 
scale problemj^0|, though these do not enable the inclu- 
sion of experimental information. Our method is efficient 
enough to enable the creation of a 647 atom model of g- 
GeSe2 using only a workstation. The method is at least 
a factor of five faster than a comparable quench from the 
melt simulation, with its inherent limitations. 

METHOD 

The obvious means to incorporate interatomic inter- 
actions into an RMC simulation is to add a constraint 
to minimize the magnitude of the force on all the atoms 
according to some energy functional or to minimize the 
total energy. For an ab initio Hamiltonian this is expen- 
sive, since Monte Carlo minimization of Eq.^ requires a 
large number of energy/force calls. Thus, we have instead 
employed a simple "self consistent" iteration scheme (in- 



dicated in Fig. (1) starting with an initial configura- 
tion Ci, minimize ^ to get C2, 2) steepest-descent quench 
C2 with an ab initio method to get C3, (3) subject the 
resulting configuration to another RMC run (minimize 
^ again), repeat steps (2) and (3) until both the MD 
relaxed model and RMC models no longer change with 
further iteration. In this paper, we limit ourselves to the 
first term in Eq. 1 (the experimental static structure 
factor), though additional constraints certainly could be 
employed. For the RMC component of the iteration, we 
make the conventional choice of using Monte Carlo for 
the minimization. This is simple and does not require 
gradients (and thus allows the use of non-analytic terms 
in Eq.nQ, if desired). 

We emphasize that our method is flexible. It's logic 
suggests that one should include whatever experimental 
information is available. In this paper we limit ourselves 
to the pair-correlation functions. In principle, other ex- 
periments could be included as well. These might be 
costly to include (for example to compel agreement with 
the vibrational density of states, the dynamical matrix 
would be required at each iteration). The method is 
equally suited to fast empirical potentials, which would 
allow studies of very large models. It is also possible to 
force a close fit to some restricted range of data, and a 
less precise fit elsewhere if desired. Our scheme also pro- 
vides insight into the topological signatures of different 
constraints (experimental or otherwise). Chemical and or 
topological constraints could also be maintained as part 
of the RMC iteration. 

Our method can be understood as a way to mini- 
mize an effective potential energy function Veff{R) = 
V{R) +AC,{R), where V{R) is the potential energy of the 
configuration (denoted by i?), A > 0, and C is a non- 
negative cost function enforcing experimental (or other) 
constraints as in Eq. ^ Empirically, we find that it is 
possible to find configurations that simultaneously ap- 
proximately minimize both terms (which implies that the 
choice of A is not very important). It is also clear that our 
method is really statistical: in general one should gener- 
ate an ensemble of conformations using ECMR. For ad- 
equately large models, self averaging can be expected; in 
this study of large (647 atom) models of g-GeSe2 we find 
similar results for two runs; for small systems a proper 
ensemble average is required. 

The method is new and as such needs to be studied and 
developed in a number of ways. Nevertheless, we show in 
this paper that it is relatively easy to model a particularly 
challenging material with significant advantages in both 
experimental plausibility of the model and computational 
efficiency of the algorithm. 
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FIG. 1; Flow diagram for the "Experimentally Constrained 
Molecular Relaxation" method of this paper. 



APPLICATION TO GLASSY GESE2 

We apply ECMR to glassy GeSe2, a classic glass 
forming material with challenging physical and tech- 
nical issues: (1) it displays nanoscale order: a "first 
sharp diffraction peak" (FSDP) is observed in neutron 
diffraction measurement, (2) the packing of GeSe tetra- 
hedra involves both edge- and corner- sharing topolo- 
gies; (3) the material has interesting photoresponse (un- 
derstanding of which requires the electronic structure), 
(4) the material is difficult to simulate with ab initio 
techniques [U [H, E [3 . The model used in our calcu- 
lation consists of 647 atoms in a cubic box of size 27.525 

A. 

In the nomenclature of Fig. ^ Ci is obtained by con- 
straining the coordination number (2 for Se, 4 for Ge) 
and the bond-angle distribution of Se-Ge-Se to an ap- 
proximate Gaussian with an average bond angle 109.5°. 
The initial network was "generic" and included none of 
the detailed local chemistry of Ge and Se aside from the 
coordination and chemical ordering (bond angles were 
not constrained in RMC loops). Equal weighting was 
used for all experimental points in this paper. Using the 
method of isotopic substitution, Salmon and Petri 
were able to separately measure the three (Ge-Se, Ge- 
Ge and Se-Se) partial structure factors of g-GeSe2. We 
jointly enforced all three partials (in real space) in the 
RMC component of the loop in Fig. ^ The MD relax- 
ation was done with FIREBALL|16|. It was found that 
after the fourth iteration, S(Q) hardly changed. In Table 
I, we show the average force per atom at the beginning 
of each call to MD relaxation; good convergence is ob- 
served. Subsequent discussion in this paper is for the 



last step of the MD, with forces less than 1 x 10~'^eV/A. 
it was not obvious to us in the beginning that RMC and 
first principles interatomic interactions could be made 
"self-consistent" , but for this system at least, reasonable 
convergence is possible. It is likely that some initial con- 
formations Ci will get "stuck" and require a new start, 
but we have not encountered difficulty with this yet. 

Structure 

In Fig. El we compare the RMC, ECMR and experi- 
mental structure factors. Here, the RMC model is that 
obtained by starting with the generic Co configuration, 
and forcing agreement on the experimental S{Q) (with- 
out any other constraints). While the agreement is very 
good, it is not perfect. This is to be expected for three 
reasons: (1) consistency between data and Hamiltonian 
is never exact; (2) our cell contains 647 atoms, which is 
compared to the thermodynamic limit and (3) we chose 
to constrain our model using real space data, which in- 
volves Fourier transforms and windowing (this introduces 
only small errors in this data set). In Fig. 2 we highlight 
the differences between experiment j 1 5j . a quench from 
the melt model0] and the new ECMR model. In the 
inset of Figure 2, we also illustrate the level of agree- 
ment using a pure RMC approach, which is similar to 
the ECMR result and notably better than quench from 
the melt. For reference, we have reproduced the full par- 
tial structure factors elsewhere 0. 

Note in Fig. |21 that the first sharp diffraction peak 
(FSDP) is well reproduced, (very close in width and cen- 
tering, and much improved from all previous models in 
height). Moreover, as for our "Decorate and Relax" (DR) 
methodp7j|. the large Q structure closely tracks experi- 
ment (unlike the experience for quench from the melt 
models which are too liquid-like and therefore decay too 
rapidly for large Q). These desirable features are of course 
"built in" ; we show here that the ECMR method does 
preserve every important feature of the structure of the 
glass manifested in S(Q). 

An important indicator of network topology and 
medium range order of GeSe2 glass is the presence 
of edge-sharing and corner-sharing tetrahedra. Raman 
spectroscopy [3| and neutron diffraction studies 
have indicated that 33% to 40% of Ge atoms are in- 
volved in edge sharing tetrahedra. The fraction in our 
model is found to be 38%. This was not "built in" to our 
modeling, and is a pleasing prediction arising from the 
procedure. We also have observed that 81% of Ge atoms 
in our model are four-fold coordinated of which approx- 
imately 75% form predominant Ge-centered structural 
motifs Ge{Sex)4 while 6% are ethane-like Ge2{Sei)e 
units. The remaining Ge atoms are three- fold coordi- 
nated and are mostly found to be bonded as Ge-Sea 
units. On the other hand, the percentage of two-, three- 
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FIG. 2: Neutron-weighted static structure factor, comparing 
ECMR model, experiment and a quench from the melt 
made with the same Hamiltonian used with ECMR,[l3j| . Inset: 
blowup of small-Q region showing initial RMC model (eg, 
enforcing experimental structure factor, but without ECMR 
iterations), experimen t 1 15l and quench from the melt model 
due to Cobb et al.\L^\l^. The first sharp diffraction peak is 
closely reproduced by ECMR and RMC, and is present but 
weak in the quenched model. 



TABLE I: The convergence of ECMR described in the text. 



ECMR iteration 



Average force/atom (eV/A) 
^■3 



2.242 X 10" 
7.365 X 10"^ 
6.518 X 10"* 
5.019 X 10"* 
4.773 X 10"'' 
4.903 X 10"'' 
4.686 X 10"" 
4.642 X 10"" 



and one-fold coordinated Se atoms are 72%, 18% and 
10% respectively. Mossbauer experiments, vi^here Sn was 
used as a Ge probe 0|, estimated that the fraction of 
Ge involved in dimers is 16% which is again in favor of 
our model. 

By integrating partial radial distribution functions 
via Fourier transform of structure factors Petri and 
Salmon [isj obtained nearest neighbor coordination num- 
bers noe-Ge = 0.25, nse-Se = 0.20, and nce-Se = 3.7 
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FIG. 3: The electronic density of states (Gaussian-broadened 
Kohn-Sham eigenvalues) for ECMR model of GeSe2, along 
with the RMC model (not using ab initio information) and a 
"decorate and relax" (DR) model made with the same Hamil- 
tonian (inset). The XPS0| and IPES 22] data show the oc- 
cupied (valence band) and unoccupied (conduction band) part 
of the spectrum. See Table-II for numerical comparison of the 
peaks. The Fermi level is at E=0. Both DR and ECMR re- 
produce the state density closely, while the RMC model lacks 
an optical gap. 



that corresponds to average coordination number n = 
2.68. The corresponding values from our model are: 
riGe-Ge = 0.17, nse~Se = 0.30, nQe-Se = 3.68, and n = 
2.66. The partial and total coordination numbers, there- 
fore, agree well with experiments (as expected) and are 
consistent with the 8-N rule which predicts n = 2.67. The 
percentage of homopolar bonds present in our model is 
found to be about 6.2 % which is again very close to the 
value 8 % noted by Petri and Salmon Ts]. 



Electronic density of states 

Having studied structural properties, we now briefly 
analyze electronic properties of our model. Since struc- 
tural and electronic properties are intimately related, an 
examination of electronic density of states provides an 
additional test of the validity of the model which is de- 
rived jointly from structural information and a suitable 
interatomic interaction. The electronic density of states 
(EDOS) is obtained by convolving each energy eigenvalue 
with suitably broadened Gaussian. The ECMR EDOS in 
the inset of Fig.|31agrees quite well with experimental re- 
sults obtained from x-ray photo-emission [21] (XPS), in- 
verse photo-emission spectroscopy '2^ (IPES) and ultra- 
violet photo-emission spectroscopy {23j (UPS) measure- 
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ments as well as with those obtained in recent theoretical 
studies [11 [Hill 113. It is remarkable that the Kohn- 
Sham eigenvalues (obtained in the Harris approximation) 
agree so well with the photoelectron spectroscopy [2^, 
particularly as the energy-dependent matrix element is 
not included in the calculation. The substantial splitting 
between the first two peaks of the valence bands named 
the Ai and A2 peaks is also well-pronounced. The po- 
sition of the principal peaks obtained from the different 
models and experiment are tabulated in Table II. The 
similarity of experiment and theory suggests the util- 
ity of a study of the Kohn-Sham eigenvectors to enable 
atomistic identification of defects and bands illustrated 
in Fig. El 

We also compare the EDOS for the RMC model. The 
RMC model does very poorly, without even showing an 
optical gap, despite the excellent static structure factor 
(obtained by construction). By contrast, our DR and the 
quench from the melt model (not shown) are very close to 
experiment and ECMR. As the coordination and chem- 
ical order is correct in the RMC model, the lack of an 
optical gap originates in an unrealistic bond angle distri- 
bution in the RMC model (something very similar hap- 
pens in a-Si if only S{q) (and no bond angle constraint) 
is used to form the modelj^. This result emphasizes the 
need to compute the density of electron states as an im- 
portant gauge of the credibility of a model. 



Vibrations 

It is useful to also examine the vibrational density of 
states (VDOS) of our ECMR model due to the close re- 
lationship to its atomic-scale structure and its dynami- 
cal properties. The VDOS was reported elsewhere 0|. 
Comparing our VDOS with experiment obtained by in- 
elastic neutron-scattering , the spectrum exhibits the 
same features with somewhat better resolution than re- 
sults we reported in Ref. Three bands can be dis- 
tinguished: a low energy acoustic band involving mainly 
extended interblock vibrations and a high energy optic 
band consisting of more localized intrablock vibrations. 
The two main bands are clearly separated by the tetra- 
hedral breathing {Ai-Aic) band. The overall agreement 
is quite reasonable, including a resolved Ai "companion" 
mode "Ale". 

In Figure 4, we compare the vibrational density of 
states of the ECMR model with experiment jl^l and for 
completeness our decorate and relax model including 
648 atoms along with the ECMR model. We do not 
present the RMC result, as the system is not at equi- 
librium according to FIREBALL, which would therefore 
lead to many imaginary frequencies in the density of 
states. While generally DR and ECMR are quite sim- 
ilar, we note some difference in the tetrahedral breathing 
Ai band (near 25 meV), including a slightly different Ai- 
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FIG. 4: Vibrational density of states computed from dynam- 
ical matrix for 648 atom models and exDeriment|27j . Nomen- 
clature similar to Fig. 3. 



TABLE II: The positions of the Ai, A2, A3 and B peaks in 
the EDOS of GeSe2 glass compared to experimental results 
0. 



(eV) 


Ai 


A2 


As 


B 


Present work 


-1.55 


-3.0 


-4.6 


-7.4 


Experiment [2^ 


-1.38 


-3.0 


-4.6 


-7.8 


Melt and quench 


-1.4 


-2.7 


-4.6 


-7.0 


Decorate and relaxfl2] 


-1.36 


-2.8 


-4.5 


-7.2 



Ale splitting. This is probably because the ratio of edge 
to corner sharing tetrahedra is different (« 29% of Ge 
atoms are involved in edge sharing tetrahedra in the DR 
model). This suggests that the VDOS has some sensitiv- 
ity to medium range order (tetrahedral packing) in this 
glass. 



CONCLUSION 

In summary, we have proposed a new method which en- 
ables the inclusion of a priori information (experimental 
or otherwise) into molecular simulation. We have shown 
that the method is effective for a challenging material, 
g-GeSe2. 
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